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Abstract 



o 

VO We present a nonlocal formulation of contact mechanics that accounts for the interplay 

^"^ of deformations due to multiple contact forces acting on a single particle. The analytical 

^^ formulation considers the effects of nonlocal mesoscopic deformations characteristic of 

confined granular systems and, therefore, removes the classical restriction of independent 
Q contacts. This is in sharp contrast to traditional contact mechanics theories, which are 

^ strictly local and assume that contacts are independent regardless the confinement of the 

^ particles. For definiteness, we restrict attention to elastic spheres in the absence of gravi- 

atational forces, adhesion or friction. Hence, a notable feature of the nonlocal formulation 
I is that, when nonlocal effects are neglected, it reduces to Hertz theory. Furthermore, we 

^T^ show that, under the preceding assumptions and up to moderate macroscopic deforma- 

C tions, the predictions of the nonlocal contact formulation are in remarkable agreement 

with detailed finite-element simulations and experimental observations, and in large dis- 
agreement with Hertz theory predictions — supporting that the assumption of independent 
contacts only holds for small deformations. The discrepancy between the extended theory 
presented in this work and Hertz theory is borne out by studying periodic homogeneous 
systems and disordered heterogeneous systems. 



<N 

> 

o 

"^ 1 Introduction 

O The understanding of confined granular systems is crucial for many fields of science and en- 

^ gineering. Compaction of granular media plays a relevant role in the design, optimization 

and control of many pharmaceutical, detergent, food, ceramic and metallurgical manufactur- 
. ^ ing processes, so much so that mechano-chemical properties of powder compacts have direct 

/\ impact on the end-product performance (see, e.g., Alderborn and Nystrom [ ], and references 

^ therein). In the context of earth sciences, the prediction of macroscopic behavior from micro- 

scopic details of densely packed systems is a central focus of soil mechanics (see, e.g., Wood 
[iv]). 

Granular systems, such as powders or soils, are made of discrete particles larger than 1- 
10 fim that interact with each other through contact forces. Despite their apparent simplicity, 
the theoretical and computational modeling of granular media remains an active area of re- 
search to this day. Since the macroscopic behavior of these systems is fundamentally encoded 
at the granular scale, one of the challenges is to develop predictive and computationally effi- 
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cient macroscopic models based solely on the bulk mechanical properties and the geometric 
configuration of the individual particles. 

The static macroscopic behavior of confined granular systems (e.g., densification and strength- 
ening) is usually modeled using either continuum mechanics or particle mechanics approaches. 
Models based on a continuum mechanics description of the granular system rely on phe- 
nomenological constitutive models such as Cam-Clay [ ] and Drucker-Prager [6] which as- 
sume that the strength of the compact is only a function of its porosity, or on micromechan- 
ical constitutive models derived by homogenization techniques (see, e.g., Cambou et al. [ ] 
and references therein). The former are computationally efficient though only predictive after 
case-by-case experimental calibrations; the latter are mathematically challenging due to the 
heterogeneous nature of deformation at the granular scale. Models based on a particle me- 
chanics approach, generally referred to as discrete element methods (DEM) [ ], improve the 
description at the particle level but their predictability relies on the contact formulations em- 
ployed. In general, at small deformations and low relative densities. Hertz theory is used to 
describe the contact behavior between two elastic particles [ ' "], and elasto-plastic contact laws 
have been proposed to account for plastic deformations (see, e.g., Vu-Qouc and Zhang [ ]). It 
bears emphasis that contact formulations are also relevant to all micromechanical constitutive 
models. Recently, a multiscale model that seamlessly bridges the particle mechanics approach 
and the continuum mechanics approach has been proposed [1 1, 20]. This quasi-discrete model 
combines a detailed description of the granular scale with the computational efficiency typical 
of finite-element discretizations of continuum models. 

DEM remains an attractive approach because of its ability to describe the complex kine- 
matics of confined granular systems (e.g., particles rearrangement). However, the need of 
predictive contact formulations at large deformations and high relative densities has led to 
new approaches where contact of individual particles is described with finite-element models. 
Then, either all particles in the system are simulated [7, 14], or force-displacement contact 
laws are extracted for DEM calculations [ ]. While both methodologies provide an accurate 
description of particles' deformations, the latter requires a case-by-case characterization of pair 
interactions, and the former restricts the analysis to a small number of particles. 

It is interesting to note that, while the deformation of elastic spheres under contact stresses 
has been extensively studied, no particular attention has been devoted to the interplay of defor- 
mations due to multiple contact forces acting on a single particle. The general assumption made 
is that contacts between particles are independent and therefore contact forces are resolved lo- 
cally. However, the deformation of an elastic sphere under contact load is not strictly local. 
This was early recognized by Tatara [ ] who studied two elastic spheres under compression 
and proposed analytical formulae that account for the effect of the reactions on the mutual sur- 
face of contact. In the context of confined granular systems, Mesarovic and Fleck [13] studied 
elasto-plastic spheres and concluded that the assumption of independent contacts only holds in 
the early stages of compaction. Recently, Harthong et al. [ ] fitted, from finite-element simu- 
lations of simple periodic lattices, contact laws that account for nonlocal-deformation effects 
during high-density compaction of elasto-plastic particles. 

The work presented in this paper is concerned with the development of a nonlocal for- 
mulation of contact mechanics that accounts for the interplay of deformations due to multiple 
contact forces acting on a single particle. The analytical formulation considers the effects 
of nonlocal mesoscopic deformations characteristic of confined granular systems and, there- 
fore, removes the classical restriction of independent contacts. This is in the spirit of [16] 
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and is in sharp contrast to traditional contact mechanics theories, which are strictly local and 
assume that contacts are independent regardless the confinement of the particles. For definite- 
ness, we restrict attention to elastic spheres in the absence of gravitational forces, adhesion or 
friction. Hence, a notable feature of the nonlocal formulation is that, when nonlocal effects 
are neglected, it reduces to Hertz theory. Furthermore, we show that, under the preceding 
assumptions and up to moderate macroscopic deformations, the predictions of the nonlocal 
contact formulation are in remarkable agreement with detailed finite-element simulations and 
experimental observations, and in large disagreement with Hertz theory — supporting that the 
assumption of independent contacts only holds for small deformations. 

The paper is organized as follows. The analytical nonlocal contact formulation is presented 
in Section 2 and is validated, with experimental observations and detailed finite-element simu- 
lations, in Section 3. In Section 4 we investigate the discrepancy between Hertz theory and the 
extended theory developed in this work. Specifically, we study different dimensional config- 
urations, different packings and different volumetric confinements of periodic homogeneous 
systems, and the consolidation and compaction of confined disordered heterogeneous systems. 
Finally, a summary and concluding remarks are collected in Section 5. 

2 Nonlocal mesoscopic deformations in compressed elastic spheres 

The classical solution to the static contact problem of elastic spheres is given by Hertz theory. 
In the absence of gravitational forces, adhesion or friction, the theory considers linear-elastic 
isotropic smooth spherical particles that deform, under normal contact forces, to accommodate 
a flat contact surface [9]. Additionally, the size of the contact surface is assumed to be much 
smaller than the dimensions and radii of curvature of the contacting bodies. Hertz theory thus 
applies to contact deformations that remain localized to a small area and that do not perturb 
other contact forces acting on the same particle. In the context of confined granular systems, 
this restriction is typically referred to as independent contacts. In this work we overcome this 
limitation by accounting for the interplay of deformations due to multiple contact forces acting 
on a single particle. Specifically, we consider the contribution to the contact interface of the 
nonlocal mesoscopic deformations induced by all other contact forces acting on the particles. 
This goal is achieved by invoking the principle of superposition and solving the boundary-value 
problem of contact stresses acting on a small region of an otherwise traction-free elastic sphere 
supported at its center. Figure 1 depicts a particle under general loading-conditions and the 
decomposition of the elastostatic problem into four single-loaded boundary-value problems. 
A schematic of the deformed configuration of an elastic sphere under the action of contact 
stresses and supported at its center is also shown in the figure — albeit under large deformations 
for illustrative purposes. 

The exact solution of the problem depicted in Figure Ic was recently reported by Zhupan- 
ska [ ]. She uses a general solution for the axisymmetric problem of an elastic sphere and 
a dual series equation approach to reduce the original boundary- value problem to a Fredholm 
integral equation of the second kind, which is then solved numerically. The analytical solution 
is achieved without replacing the contact area of the sphere by an elastic half-space (relaxing 
Hertz's assumption of small contact area) and it allows for determining the contact pressure 
and the displacement of the traction-free surface. The results presented in [21] show that the 
spherically-distributed contact pressure predicted by Hertz theory remains accurate for rela- 
tively large contact areas. The aim of our work, however, is to solve an approximate problem 
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Figure 1: a) Particle within a confined granular system under general loading-conditions, b) 
Decomposition, by invoking the principle of superposition, of an elastostatic problem into four 
single-loaded boundary-value problems, c) Schematic of the deformed configuration of an 
elastic sphere under the action of contact stresses and supported at its center — albeit under 
large deformations for illustrative purposes. 

that not only allows for accounting the interplay of deformations due to multiple contact forces 
acting on a single particle, but also that admits an efficient closed-form solution. 

Next, in Section 2.1, we propose an approximate problem to the single-loaded boundary- 
value problem that is amenable to an efficient closed-form solution, and we assess its accuracy 
with respect to the exact solution reported in [" ' ]. In Section 2.2 we present a nonlocal contact 
formulation that employs this approximate solution and, therefore, by invoking the principle 
of superposition, it accounts for the interplay of deformations due to multiple contacts acting 
on a single particle. Finally, in Section 2.3, we provide some details about the computational 
implementation of the nonlocal formulation. 



2.1 Single-loaded boundary- value problem 

We consider an elastic sphere of radius R that deforms to accommodate a flat circular contact 
surface of radius a under the action of a distributed contact pressure with maximum value pm 
and supported at its center by a concentrated force P. We adopt cylindrical coordinates (2:, r), 
with the z axis aligned with the loading direction, and we denote a surface point by its angle 
9, as depicted on the left side of Figure 2. The main approximations that allow for an efficient 
closed-form solution are then twofold: 

1. Displacements and tractions of the contact surface, i.e., < 9 < arcsin(a/i?), are 
described by Hertz theory and therefore approximated by the solution of a semi-infinite 
elastic body under the action of a spherically-distributed contact pressure with maximum 
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value pm = 3P/27ra^. Thus, according to the theory of elasticity, in particular, the 
Boussinesq solution [9, 17], the compressive surface displacement Uz within the loaded 
region is given by 



uz{e) 



3P(1 



Sa^E 



{2a^ -R^sm^ie)) 



(1) 



where 9 G [0; arcsin(a/i?)], E is the Young modulus, and v is the Poisson's ratio. 



Displacements of the traction free surface, i.e., arcsin(a/i?) < ^ < tt, are approximated 
by the solution of a spherical particle under a concentrated force P — applied at the origin 
of coordinates and aligned with the z axis — and supported by a small surface traction 

qi^) = [qzAr) given by 



m = 



-3sin(e/2), 



l-2i/ 



87ri?2 \^ "^"^'"'"'' 2tan(6'/2)(l + sin(6'/2)) 
Thus, the normal displacement of a surface point is 

l + v P [-2(1 -v)- 2(1 - 2v) sin((9/2) + (7 - 



3cos(6'/2)\ 
2 J 



un{e) 



w) sm' 



\e/2)] 



AttRE 



sin(0/2) 



(2) 



(3) 



where 9 G (arcsin(a/i?) ; n] and the same displacement reference frame employed above 
is chosen (i.e., for this case, Uz{7r) -^ as R ^ oc). It is worth noting that q scales as 
ll^ll/Pm = O{o? /R^), that is the surface traction is much smaller than contact forces. 
We present in the Appendix a detailed derivation of this result together with other possi- 
ble approximate problems amenable to closed-form solutions. 
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Figure 2: Elastic sphere under the action of contact stresses and supported at its center, 
for a/R = 0.233 and v = 0.50. Left: Loading configurations of exact and approximate 
problems — the surface traction scales as ||g||/pm = O{o? /R^). Right: Surface normal dis- 
placement given by the exact solution (symbols), the approximate solution (solid curve), and 
by traditional contact mechanics theories (dashed curve). The exact solution is taken from 
Zhupanska [21]. 

Under the preceding assumptions, the displacement of surface points outside the contact 
area can be readily computed. We then proceed to compare this approximate solution with the 
exact solution presented by Zhupanska [ ]. The left side of Figure 2 illustrates the loading 
configuration of the exact boundary- value problem and the approximate problem for a/R = 
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0.233 and z/ = 0.50. It is evident that the surface traction is much smaller than contact forces 
even for a non-small contact area. Since an incompressible material, i.e., z/ = 0.50, experiences 
the largest displacements, this is the most interesting case for assessing the accuracy of the 
approximation. The right side of the figure shows the surface normal displacement given by 
the exact and approximate solutions, and by traditional contact mechanics theories — a^ = 
3PR{1 — u'^)/4:E is used in the representation [9]. The results are in very good agreement and 
put into evidence that the interplay of deformations due to multiple contact forces acting on a 
single particle is not necessarily negligible. 

2.2 Nonlocal contact formulation 

We now consider the displacement of the contact surface between two elastic spheres of radius 
Ri and R2, and material properties Ei, ui and E2, 1^2, respectively. The spheres are com- 
pressed by a general loading configuration of concentrated forces that are statically equivalent 
to an effective force P aligned with the z axis, as depicted in Figure 3. Following Hertz theory, 
the effective force is assumed spherically-distributed within a flat circular contact surface of 
radius a. Thus, a point within the contact surface verifies the following compatibility equation 



7 = ^1 



Rl 



r2 + i?2 



i?2 



r^ + wi{r) + W2{r) 



k=l,2 i=l 



i{r) (4) 



where 7 is the relative displacement along the z axis of the center of masses, and wi, W2 are 
the displacements due to the local contact pressure. By invoking the principle of superposi- 
tion, Wpi are the contributions of the nonlocal mesoscopic-deformations induced by all other 

k 

concentrated forces P^ acting on each spherical particle k. 




Figure 3: Two spheres in contact under a general configuration of multiple concentrated forces. 
Left: Schematic of the loading configuration with all forces rotated about the z axis and repre- 
sented on the zr plane. Middle: Contact loading-conditions for Hertz theory. Right: Contact 
loading-conditions for the nonlocal contact formulation proposed in this work. 

The compatibility equation (4) then provides the equilibrium configuration 7 and the con- 
tact surface radius a for a given loading configuration P^, 0^. The evaluation of Wk{r) is 
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approximated by the solution of a semi-infinite elastic body under a spherically-distributed 
contact pressure — i.e., w{t) = iX2(arcsin(r/i?)) from Equation (1). In addition, Wpi{r) is 
approximated by a single value equal to the normal displacement induced along the z axis by 
the force P^, that is Equation (3), 



l + uj, PI [-2(1 - z/fe) - 2(1 - 2z/fe) sin(g^/2 ) + (7 - Suj,) sin2(0^/2)] 
^ ^TTRkEk sin(0^/2) 



^ ____^ ^^_^ — ._. — ^ — r ..;.:: ' ^ — - — ^-^^ (5) 



where O], is the angle between the concentrated force and the z axis (see Figure 3). Finally, the 
profile of the undeformed spherical cap is replaced by the first term of its series expansion 



Ri - JrI -r^ = Ri 



1 r2 „ / a4 



2Ry^\Ri 



1 



„2 



2Ri 



r' (6) 



After introducing the preceding approximations, equilibrium configuration for 7 and a, 
given the loading configuration statically equivalent to P, can be readily solved 

/ px2/3 

7 = — - INL (7) 



El E2 J \Ri i?2 

where uh corresponds to Hertz theory and ^nl accounts for the nonlocal loading terms (i.e., 
all Wpi defined by (5)) 

k 



^k 



riH 4.\ El E2 J \Ri R2 

N,, (9) 

k=l,2 i=l 



Therefore, if nonlocal terms are neglected, the formulation reduces to Hertz theory, i.e., ^(7) = 
riHj^^^, and the elastic contact depends on the relative position of the center of masses regard- 
less the loading configuration (see the middle part of Figure 3). In sharp contrast, if nonlocal 
terms are considered, the formulation departs from Hertz theory and the elastic contact now 
depends on the loading configuration of each sphere (see the right part of Figure 3). 

The formulation presented above, i.e.. Equations (7), (8) and (9), is explicit provided the 
network of contact forces and loads are known. However, this is not the focus of this study. 
Instead, we seek for the equilibrium configuration of a system of particles given loads and 
static boundary conditions. That is, we aim at determining the relative displacements 7 that 
generate a network of forces P in static equilibrium with the loads and boundary conditions. 
Hence, we first rewrite Equation (7) as follows 

^(7) = ^if (7 + INL^^ (10) 

to obtain an explicit formulation in terms of relative displacements. Then, the solution of the 
non-linear system can be efficiently obtained by a few fixed point iterations of P, starting with 
INL = 0. It bears emphasis that jnl is a function of the contact forces and, in the context of a 
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Figure 4: Two examples that illustrate nonlocal mesoscopic-deformation effects in contact 
pair-interactions. Left: Two elastic spheres in compression behave stiffer than Hertz theory 
predictions. Right: Two elastic spheres are brought into contact by transverse loading. 



numerical iterative procedure, it has to be approximated by the values of its previous numerical 
guess. 

It is interesting to verify that Equations (10) and (9) simplify to the results obtained by 
Tatara [ ] when two elastic spheres are compressed by forces F aligned with the z axis (see 
left side of Figure 4). Then, for a loading configuration given by P/ = Pj^ = F and 9\ = 



92 = TT, the contribution of the nonlocal forces is 



INL 



F{l + Vk){'i-2vk) 



k=ia 



^TiRkEk 



(11) 



and three consecutive iterations of the force F are 



F = riif 7' 
F 



3/2 



+ 0(7' 



HHl^^'^ + 



1^7^ + Oil'") 



(12) 



F 



riHl 



3/2 



3n^ _2 



2nNL 



7 + 



2ln 



8n%^ 



H ^5/2 + Oil' 



with riNL = F/lNL' Since the first iteration corresponds to jnl = 0, the value of F coincides 
with Hertz theory prediction. These are indeed the formulae obtained by Tatara in his work^ . 

It is also interesting to observe that there are loading configurations that result in compres- 
sive forces P > with no relative displacement of the particles, i.e., 7 = 0. This is the case of 
the example depicted on the right side of Figure 4 where the loading configuration is given by 
Pi = Pf = Pj ^P^^F and e\^el^e\^el^ 7r/2. Specifically, the reaction P at the 
supports that enforce 7 = is given by 



riH 

3/2 



i?3/2 



where unl = F/^nl is 



1 

riNL 



E 

k=l,2 



(1 + Vk){^\f2 - 4 - 4:V2vk + St^fc) 
Ai^RkEk 



(13) 



(14) 



It is worth noting that Hertz theory predicts P = for this configuration, regardless the value 
ofF. 



'Besides a (perhaps typographic) error in the third term of the series expansion. 
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2.3 Implementation of the nonlocal contact formulation 

The implementation of the nonlocal contact formulation presented in this section requires the 
solution of a set of non-linear equations. The set of non-linear equations for a heterogeneous 
granular system is uniquely defined by the material properties and geometry of the individual 
particles, a simplicial complex X comprised by a set of vertices Eo{X) that represent the 
center of mass of each particle and by a set of edges Ei{X) that represent plausible contact 
pair-interactions, and the coordinates qe^ of all vertices Ca G Eo{X). The left side of Figure 5 
depicts a simplicial complex and the underlying heterogeneous granular system; additionally, 
vertices a and b, and edge ab are highlighted. 




.• 



^NL^lNL+<i 



Figure 5: Left: Heterogeneous granular system and its corresponding simplicial complex. Ver- 
tices represent centers of mass and edges represent plausible contact pair-interactions (edge ab 
is highlighted). Right: Incident edges to edge ab are grouped by those connected to vertex a 
and those connected to vertex b, that is those contact pair-interactions acting on particle a and 
those acting on particle b. 

The solution of the non-linear system can be efficiently obtained by a few fixed point 
iterations of the nonlocal terms {7^7^} defined by equation (9). An iteration reduces to two 
nested-loops: (i) compute from equation (10) all contact forces P^"^ associated with edges Ca, 
(ii) update nonlocal terms 7^^ with the contributions We^ induced in all incident contact pair- 
interactions 6/3 (the right side of Figure 5 illustrates all incident edges to edge ab). Finally, 
contact forces are determined from converged nonlocal contributions {7^^} and assembled 
into global forces F (i.e., F is a vector that collects the out-of-balanced force of each particle 
in the system). Thus, the nonlocal contact formulation only adds a computational overhead 
cost of complexity OinNc), where Nc is the number of contact pair-interactions and n is the 
average number of neighbouring particles. The implementation is summarized in Algorithm 1. 
It is worth noting that this is the only modification required for implementing the nonlocal 
contact formulation in existing codes. 

3 Validation of the nonlocal contact formulation 



The validation of the nonlocal contact formulation presented in this work is twofold. First, 
we compare the analytical predictions with experimental results from Tatara [K] for a rubber 
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Algorithm 1 Nonlocal contact formulation. Computation of global forces. 



Require: {7^°^ }, TOL, and the simplicial complex X. 



9 

10: 
11; 
12: 
13; 

14; 
15; 
16; 

17; 
18; 
19; 
20; 
21; 
22; 
23; 
24; 
25; 
26; 
27; 
28; 
29; 



/* Compute nonlocal contributions {^y^j^} */ 
Error ^ TOL 

{'1n\} ^ {0} 
while Error > TOL do 

/* Loop over all contact pair-interactions e^ */ 
forec, ^Ei{X)Ao 

if 7^- +^7^"^ >Othen 

/* Compute contact force from (10) */ 

P^- ^n^(7^"+S^"^)^/^ 

/* L6X9/7 6>v^r all contact pair-interactions e^ G £'i(X) incident to Ca */ 

for 6/3 G Incident (ec,) do 

/* Update 7^^ w///z //z^ nonlocal displacement 

induced by contact force P^"^, i.e., We^ from (5) */ 

end for 
end if 
end for 

/* Compute a measure of convergence */ 

Error ^ H^+^tvl - ^7nl\\ 

k^ k + 1 
end while 

/* Assemble global force F using converged nonlocal terms {^7^7^} */ 
forec G £;i(X)do 

if 7^- +^7nl >Othen 

P^- ^n^(7^"+S^"^)^/^ 

F ^ ASSEMBLE(ec,, P"^"^) 

end if 
end for 
return F 
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sphere pressed between two rigid plates. Second, we perform detailed finite-element simula- 
tions for different loading configurations and compare contact forces with predictions of the 
nonlocal formulation. 

The experimental observations reported by Tatara [ ] correspond to rubber spheres of 
radius R = 0.01 m and elastic properties E = 1.85 MPa, u = 0.46. The spheres were pressed 
between two rigid plates in a tensile-compression test machine and, by exploring deformations 
up to 40% of the undeformed diameter, the author measured applied loads with a load-cell. 
Thus, applied load versus deformation curves were recorded. Tatara also reported that neither 
hysteresis nor permanent deformations were observed during the experiments. 

A Lagrangian finite-element formulation is adopted and finite deformations are considered 
p^]. Due to geometric and loading symmetries only one-eighth of the sphere is modeled. The 
mesh comprises 324,000 10-node tetrahedral isoparametric elements and 461,404 nodes. The 
three-dimensional solid is characterized by a strain-energy density of the form 



W{F) 



A-^M 



{\ogjf-fi\ogJ+^tr{\/^^V^) 



(15) 



which describes a neo-Hookean material extended to the compressible range. In this expres- 
sion, (f is the deformation mapping, V(f is the deformation gradient tensor, J = det(V(/^), and 
A and /x are the Lame constants — the values employed in the simulations correspond to those 
reported by Tatara [ ]. Rigid frictionless boundary constraints are enforced by a penalization 
method, whereas symmetry conditions are imposed by restraining the corresponding degrees 
of freedom. The numerical solution is presumed ostensibly converged with respect to the mesh 
size — meshes with higher levels of refinement were used by way of reference. By way of 
illustration. Figure 6 shows a coarse mesh of one-eighth of a sphere. 





Figure 6: Finite-element mesh of one-eighth of a sphere — a coarse mesh comprised by 40,500 
tetrahedral elements and 61,504 nodes is depicted. Left: Initial mesh. Right: Deformed mesh 
compressed in one direction and constrained in two other directions. 



Figure 7 summarizes the results of the validation in an applied-load- versus-deformation 
plot. Hertz theory predictions depart from experimental observations at very small deforma- 
tions. In sharp contrast, the detailed finite-element simulation and the analytical predictions 
of the nonlocal contact formulation presented in this work are in excellent agreement with ex- 
perimental measurements over the whole range of deformations. It bears emphasis that the 
computational effort required to evaluate a solution with the nonlocal formulation, or with 
Hertz theory, is negligible in comparison with the effort of solving the finite-element model. 
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Figure?: Validation of the nonlocal contact formulation. Left: Schematic of the loading config- 
uration. Right: Applied load versus deformation curves. Hertz theory prediction (black-dashed 
curve), nonlocal contact formulation results (black curve), finite element solution (grey curve), 
and experimental measurements from Tatara (1989) (five-pointed stars). 

In order to asses the accuracy of the nonlocal contact formulation in configurations where 
nonlocal mesoscopic deformations play an even larger role, we study a spherical elastic particle 
compressed in one direction and constrained in other directions. Due to the lack of experimen- 
tal observations, we perform detailed finite-element simulations and compare contact forces 
with predictions of the nonlocal formulation. Specifically, Figure 8 shows the applied load 
and the lateral reactions versus deformation of configurations with one set of lateral constraints 
and with two orthogonal sets of lateral constraints. The lateral reaction is a consequence of 
the lateral geometric confinement of the systems. Hertz theory entirely disregards such effect 
and predicts no lateral reaction forces. The present formulation, in contrast, is designed to 
account for nonlocal mesoscopic-deformation effects and therefore lateral reactions naturally 
result from the compressive loading conditions. The agreement between predictions of the 
finite-element model and the nonlocal contact formulation is again remarkable. 
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Figure 8: Validation of the nonlocal contact formulation. Left: Schematic of a spherical par- 
ticle compressed in one direction and constrained in other direction. Right: Applied load and 
lateral reaction versus deformation curves for a spherical particle compressed in one direction 
and constrained in one/two directions. Hertz theory prediction (black-dashed curve), nonlo- 
cal contact formulation results (black curve), and finite element solution (grey curve). Both 
loading configurations are intentionally depicted with the same color scheme. 
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It bears emphasis that despite the good agreement at large effective deformations, the exten- 
sion of the nonlocal contact formulation to nonlinear material models and to finite deformations 
remains a topic of interest. 

4 Macroscopic behavior of confined granular systems 

4.1 Granular crystals 

Granular crystals, or highly packed granular lattices, allow for a systematic understanding of 
nonlocal mesoscopic-deformation effects in the macroscopic behavior of granular systems. 
Specifically, we study different dimensional configurations, different packings and different 
volumetric confinements of periodic homogeneous systems, i.e., all particles are of same size 
and material properties. By way of example, and for later reference, we illustrate in Figure 9 
a one-dimensional chain, a bi-dimensional square lattice, and a three-dimensional cubic lattice 
of spherical beads. It is worth noting that the loading configurations discussed in Section 3 are 
equivalent to the granular crystals depicted in Figure 9. The experimental setup is equivalent to 
the one-dimensional configuration, and the configurations constrained in one and two orthogo- 
nal directions correspond to bi-dimensional and three-dimensional arrangements, respectively. 



QOQQ( 






a) ID configuration b) 2D configuration c) 3D configuration 

Figure 9: Different granular crystals assemblies. The bi-dimensional and three-dimensional 
configurations are periodic in all directions orthogonal to the load. 

The discrepancy between Hertz theory and the extended theory presented in this work 
evidently depends on the granular assembly and on the loading conditions. Hence, we restrict 
attention to deformations smaller than 10% and we quantify the discrepancy between the two 
formulations by 

where na^^l^ is the prediction of Hertz theory and ^(7) corresponds to the result of the 
nonlocal theory, i.e., the solution of Equations (10) and (9). Then, for a one-dimensional 
configuration. Hertz theory differs from the nonlocal contact formulation by 
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where e = 2^. For bi-dimensional and three-dimensional configurations, the biaxial and triax- 
ial loading conditions introduce a new term of order 0{e^/^) and therefore the first two terms 
in the error expansion are unperturbed, i.e.. 



^2D,3D 



1 

2^ 



3-2z/ 



1 



jy 



eV2 + 



247r2 V 1 



3-2i/ 



u 



e + 0(6^/4) 



(18) 



It bears emphasis that the discrepancy between the nonlocal contact formulation and Hertz 
theory depends strongly on the deformation e and the Poisson's ratio u, regardless the value of 
the Young modulus and the size of the particles. It also depends weakly on the dimensional 
configuration of the granular crystal. This last statement is in agreement with the numerical 
results presented in Figure 8. The error bounds presented above are illustrated in Figure 10 for 

ue (0,0.50). 
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Figure 10: Discrepancy between the nonlocal contact formulation and Hertz theory for 
u G (0,0.50). Left: One-dimensional configuration (Figure 9a). Right: Three-dimensional 
configuration (Figure 9c). 

The weak dependency on the dimensional configuration of the granular crystal suggests 
the study of square and cubic lattices loaded in different directions, e.g., the configuration 
depicted in Figure 11, and of different assemblies, e.g., the close-packed configuration depicted 
in Figure 12. Similarly, it is interesting to investigate the role of initial confinement on the 
macroscopic behavior of the system, e.g., the preloaded configuration depicted in Figure 13. 
These three cases are examined next in turn. For a displacement 7 in the direction of loading, 

and e = 



V2R 



deformations e of rotated and close-packed configurations are defined by e — ^^ 
^, respectively. 

We note from Figures 10 and 11 that nonlocal mesoscopic-deformations in lattices of in- 
compressible materials (i.e., u -^ 1/2) are more relevant for rotated configurations. However, 
highly compressible materials (i.e., u ^ 0) are less sensible to nonlocal effects in rotated 
square lattice assemblies. Granular crystals comprised by rotated square or cubic lattices are 
difficult to realize experimentally and therefore close-packed configurations become a more 
attractive approach. Figure 12 shows that the discrepancy between the nonlocal contact for- 
mulation and Hertz theory is again larger for incompressible materials and smaller for highly 
compressible particles, when compared with the previous configurations. It is also interesting 
to observe that the nonlocal formulation predicts the formation of a gap between horizontal 
neighbors of close-packed configurations of highly compressible beads. 
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Figure 11: Left: Bi-dimensional granular crystal rotated by 7r/4 — the configuration is peri- 
odic in the direction orthogonal to the load. Right: Discrepancy between the nonlocal contact 
formulation and Hertz theory for v G (0, 0.50). 



iii.iiii..liiiii 



, 30 



o 

CO 20 

Q. 

CD 

o 15 




CO 



10- 



tttnttntti! 









\y^ 








^y^ 






^ 


^ 




v=0.50 


/^ 






y^ 






y 


^ 




^ 


y 


^ 


^-^^ 


v=0.00 


C^-'-'^ 









4 6 

Deformation f%l 



Figure 12: Left: Close-packed bi-dimensional granular crystal — the configuration is periodic 
in the direction orthogonal to the load. Right: Discrepancy between the nonlocal contact for- 
mulation and Hertz theory for v G (0, 0.50). 
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Figure 13: Left: Bi-dimensional granular crystal preloaded by eg in both directions. Right: 
Discrepancy between the nonlocal contact formulation and Hertz theory for 6o = 4% and 

z/G (0,0.50). 
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The effect of volumetric confinement is shown in Figure 13. In order to confine a bi- 
dimensional granular crystal by eo in both directions, an initial load has to be applied to the 
system. Therefore, if the confined configuration is taken as the reference configuration, dis- 
crepancies between the nonlocal contact formulation and Hertz theory will be significant even 
for e = 0. Moreover, and in contrast to all previous configurations, it is evident in the fig- 
ure that the error is not order (9(e^/^) at small deformations. This different behavior suggests 
an additional validation of the nonlocal formulation with a detailed finite-element simulation. 
Figure 14 shows the applied load, and lateral reaction, versus deformation curves for a bi- 
dimensional system with an initial confinement of eo = 4% and subsequent compression in 
one direction. The agreement between the detailed finite-element simulation and the analytical 
predictions of the nonlocal contact formulation is remarkable. 




4 6 
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Figure 14: Applied load, and lateral reaction, versus deformation for a bi-dimensional granular 
crystal preloaded by eg = 4% configuration. Hertz theory prediction (black-dashed curve), 
nonlocal contact formulation results (black curve), and finite element solution (grey curve). 
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Figure 15: Body centered cubic granular crystal under hydrostatic compaction and die com- 
paction. 



Finally, we study the compaction of a three-dimensional configuration. Specifically, the 
behavior of a body centered cubic granular crystal under hydrostatic and die compaction is 
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investigated (see Figure 15). The relative density is adopted as measure of deformation, that is 



7T 



PR = 



V3 



[l-e,){l-ey){l-e. 



(19) 



where Cx, Cy, e^ are the deformations of the unit cell in each coordinate direction (e.g., for 
hydrostatic compaction e^^ = e^ = e^, and for die compaction e^^ = e^ = with e^ ^ 
0). Then, the undeformed system has pR = 0.68. It is worth noting that second contacts 
are activated during the compaction of a body centered cubic configuration — indeed, for die 
compaction, third contacts are also activated at large relative densities. This phenomenon 
results in an increasingly stiffer system and therefore large loads are required to achieve high- 
density compaction. Evidently, the accurate prediction of the onset of second and third contacts 
relies on the ability of the model to account for nonlocal mesoscopic-deformation effects. 

The right side of Figure 16 shows the discrepancy between the nonlocal contact formu- 
lation and Hertz theory in the prediction of the loading force required to achieve hydrostatic 
compaction of the granular system. For highly incompressible materials (i.e., u -^ 1/2), it is 
interesting to observe that the discrepancy first changes its curvature and then rapidly grows for 
relative densities larger than 0.95. This result is explained by the formation of second contacts. 
The nonlocal formulation predicts the formation of second contacts at pR = 0.95 for u = 0.45, 
whereas Hertz theory at pR = 1.04 for all Poisson's ratios. The left side of the figure shows 
the applied load versus relative density obtained with both formulations for a body centered 
cubic packing of rubber spheres of radius i? = 0.01 m and elastic properties £^ = 1.85 MPa, 
u = 0.46. The discrepancy between the nonlocal formulation and Hertz theory is evident in 
the figure. 
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Figure 16: Body centered cubic granular crystal under hydrostatic compaction. Left: Applied 
load versus relative density predicted by Hertz theory (black-dashed curve) and by the non- 
local contact formulation (black curve) — R = 0.01 m, E = 1.85 MPa, u = 0.46. Right: 
Discrepancy between the nonlocal contact formulation and Hertz theory for u G (0, 0.44). 



The right side of Figure 17 shows the discrepancy between the nonlocal contact formulation 
and Hertz theory in the prediction of the loading force required to achieve die compaction of 
the granular system. It is worth noting the presence of a corner point for values of relative 
density in the vicinity of 0.78. Again, this behavior is explained by the formation of second 
contacts. The nonlocal formulation predicts the formation of second contacts at pR = 0.780 
for u = 0.45, whereas Hertz theory at pR = 0.786 for all Poisson's ratios. The left side of 
the figure shows the applied load, and lateral reaction, versus relative density obtained with 
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both formulations for a body centered cubic packing of rubber spheres of radius i? = 0.01 m 
and elastic properties £^ = 1.85 MPa, u = 0.46. In this example, the nonlocal formulation 
predicts the formation of third contacts at pR = 0.96. As it was for the previous example, the 
discrepancy between both formulations is evident in the figure. 
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Figure 17: Body centered cubic granular crystal under die compaction. Left: Applied load, 
and lateral reaction, versus relative density predicted by Hertz theory (black-dashed curve) 
and by the nonlocal contact formulation (black curve) — R = 0.01 m, £^ = 1.85 MPa, u = 
0.46. Right: Discrepancy between the nonlocal contact formulation and Hertz theory for u G 

(0,0.44). 



4.2 Disordered granular media 

The understanding of nonlocal mesoscopic-deformation effects in the macroscopic behavior 
of disordered granular media is mostly application-driven. In this work, we focus on the con- 
solidation and compaction of confined granular systems (see, e.g., Redanz and Fleck [ ]). 
Specifically, we study the compaction of a powder bed comprised of 125 spherical particles 
of two different materials. The particle size and elastic properties of the two powders are: 
i) El = 19 GPa, ui = 0.30, a log-normal distribution of particle sizes with cut-offs at ra- 
diuses 80 pm and 75 /xm (typical properties of, e.g., cellulose), ii) E2 = 7 GPa, z/2 = 0.20, a 
log-normal distribution of particle sizes with cut-offs at radiuses 167 ^m and 125 fim (typical 
properties of, e.g., lactose). The die has a square cross-section of 1 mm x 1 mm. The initial 
condition of the un-compacted granular bed is 1 mm in height and it is numerically generated 
by means of a ballistic deposition technique [ 0]. Thus, the initial granular bed has a relative 
density of pR = 0.42. The left side of Figure 18 shows the un-compacted granular bed, mate- 
rials are depicted with different shades of grey. The right side of the figure shows a snapshot 
of the powder bed at pi^ = 0.56 with some particle rearrangement typical of the consolidation 
process. 

The left side of Figure 19 shows the applied load, and lateral reactions, versus relative 
density predicted by Hertz theory and by the nonlocal contact formulation. It is interesting 
to note that both formulations predict a similar deformation pattern during the consolidation 
process, that is relative density is increased by particle rearrangement as opposed to particle 
deformation. However, after consolidation is completed, mechanical compaction is required 
to increase relative density. Therefore, mesoscopic-deformation effects start to play a relevant 
role in the behavior of the confined granular system and the predictions of Hertz theory depart 
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Figure 18: Left: Initial granular bed with a relative density of pR = 0.42. Materials are 
depicted with different shades of grey. Right: Granular bed during consolidation, pR = 0.56. 
Particle rearrangement is evident. 
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Figure 19: Consolidation and compaction of the disordered heterogeneous granular bed de- 
picted in Figure 18. Left: Applied load, and lateral reaction, versus relative density predicted 
by Hertz theory (black-dashed curve) and by the nonlocal contact formulation (black curve). 
Right: Granular bed for pR = 0.80. 
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from those of the nonlocal contact formulation. Again, as it was the case for highly packed 
granular lattices, the discrepancy between both formulations is significant and evident in the 
figure. Finally, the granular bed for pji = 0.80 is depicted on the right side of Figure 19. 

5 Summary and discussion 

We have developed a new nonlocal formulation of contact mechanics that accounts for the in- 
terplay of deformations due to multiple contact forces acting on a single particle. This property 
is achieved by considering the contribution to each contact interface of the nonlocal mesoscopic 
deformations induced by all other contact forces acting on the particles. As a consequence, the 
formulation presented here overcomes the unrealistic assumption that contacts are indepen- 
dent regardless the confinement of the granular system — typical of standard contact mechanics 
theories. For definiteness, we have restricted attention to elastic spheres in the absence of grav- 
itational forces, adhesion or friction. Hence, a notable feature of the nonlocal formulation is 
that, when nonlocal effects are neglected, it reduces to Hertz theory. We have shown that, up 
to moderate macroscopic deformations, the predictions of the nonlocal contact formulation are 
in remarkable agreement with detailed finite-element simulations and experimental observa- 
tions, and in large disagreement with Hertz theory. The discrepancy between Hertz theory and 
the extended theory presented in this work has been borne out by studying selected confined 
granular systems. Specifically, we have investigated different dimensional configurations, dif- 
ferent packings and different volumetric confinements of periodic homogeneous systems, and 
the consolidation and compaction of confined disordered heterogeneous systems. For all stud- 
ied cases, the discrepancy between both formulations suggests that mesoscopic-deformation 
effects play a relevant role in the behavior of the granular system. Characteristically, for nearly 
incompressible materials (i.e., v -^ 1/2), discrepancies of 40% are observed in periodic bi- 
dimensional granular configurations under 10% of macroscopic deformations, discrepancies 
of 200% are predicted for nearly full compaction of body centered cubic granular crystals. 
Discrepancies in the order of 50% are observed during consolidation and compaction of an 
heterogeneous and disordered granular bed comprised by relatively compressible materials 
(i.e., V ^ 1/4). 

We close by pointing out some limitations of our approach and possible avenues for exten- 
sions of the formulation. 

First, we have restricted attention to linear-elastic isotropic smooth spherical particles. The 
systematic extension of the nonlocal contact formulation to nonlinear material models, such as 
neo-Hookean and elasto-plastic constitutive relations, to finite deformations, and to bonding 
interactions are worthwhile directions of future research. 

Second, it is clear that the presented approximations to the single- valued boundary value 
problem are not the only — perhaps even the best — that are amenable to an efficient closed-form 
solution. The investigation of other approximation strategies and the determination of the most 
accurate approximate solution are promising research directions. 

Third, although the application of the nonlocal contact formulation to dynamic problems 
is straightforward, we have only focused on static analyses. The effect of nonlocal mesoscopic 
deformations in the dynamic macroscopic behavior of confined granular systems remains as an 
interesting topic to study. 

Fourth, it is worth to note the similarity with which discrete granular and atomistic sys- 
tems are theoretically described despite their different time and length scales. Specifically, 
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iteratomic forces derived from pair potentials in discrete atomistic simulations are the counter- 
parts of the contact forces obtained from traditional contact mechanics theories. In the context 
of our work, the nonlocal formulation has conceptual similarities with embedded-atom poten- 
tials which are composed by pairwise terms and many-body terms [ ]. Indeed, the many-body 
terms account for the superposition of electron densities of the surrounding atoms. However, a 
rigorous analysis of such correlation is desirable, if beyond the scope of this paper. 
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Appendix 

We proceed to present two approximate problems to the boundary-value problem depicted in 
Figure Ic that admit efficient closed-form solutions. 





a) Exact boundary-value problem. 



b) Case I. 



c) Case 11. 



Figure 20: Elastic sphere under the action of contact stresses and supported at its center, for 
a/R = 0.233 and u = 0.50. Loading configuration of the exact boundary-value problem (a), 
and the two proposed approximate problems referred to as Case I (b) and Case II (c). The 
surface traction scales as ||g||/pm = O{o? /R^). 



Case I 

The first approximate problem is the one employed in this paper and corresponds to an elastic 
spherical particle under a concentrated force P — applied at the origin of cylindrical coordi- 
nates (z, r) and aligned with the z axis — and supported by a small surface traction q given by 
Equation (2). Figures 20a-b illustrate the approximate problem and the exact problem for an 
incompressible material and a/R = 0.233, where a is the radius of the contact surface and R 
is the radius of the elastic sphere. 

In order to render the approximate problem exactly solvable, we adopt the displacement 
and stress fields of a semi-infinite elastic body B under a concentrated force P applied at the 
origin of cylindrical coordinates, i.e., the Boussinesq solution [9, 17]. We then restrict attention 
to a spherical subbody >Bsph(i?) of radius R and tangent to the z-plane at the application point 
of the force (see Figure 20b). The surface points that describe the sphere, i.e., points (z, r) G 
^^Sph(i?)' ^^^ given by (2i?sin^(0/2), i?sin(0)), and the unit normal to 5>Bsph(i?) is n = 
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(— cos(^), sin(^)). Next, we replace the part of the body outside the sphere, i.e., B — l3sph{R)^ 
by the equilibrium surface traction q acting on 5>Bsph(i?), that is 
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(20) 



where the stress components correspond to the Boussinesq solution for a semi-infinite elastic 
body under a concentrated force. 

Based on Boussinesq's solution, the displacement of a surface point in 5>Bsph(i?) is given 
by 
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where E is the Young modulus and v is the Poissons ratio. The components of the stress 
field — relevant to our analysis — at such surface point are 
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Thus, after some trite manipulations, the surface traction q simplifies to 
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and the displacement of a surface point in the normal direction is given by 
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Naturally, the above equations correspond to Equations (2) and (3) presented in Section 2. 
Finally, it bears emphasis that the surface traction scales as 
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where pm is the maximum value of the contact pressure (e.g., pm = 3P/2Tra? for a spherically- 
distributed contact pressure). It is also important to note that the displacement reference frame 
employed in the derivation is the one of Boussinesq's solution (i.e., Un(7r) ^ as i? — > oo). 

Case II 

The second approximate problem aims at a more accurate representation of the contact pres- 
sure. We assume a spherically-distributed contact pressure p{r) = Pm\/l — r2/a2 applied on 



M. Gonzalez et al. 



A nonlocal contact formulation 



23 



a flat circular surface, with radius a, of an otherwise spherical elastic particle of radius R. Sim- 
ilarly to the previous case, the sphere is supported by a small surface traction q that scales as 
ll^ll/Pm = O{o? /R^). Figure 20c illustrates the approximate problem for an incompressible 
material and a/R = 0.233. 

Following the strategy used in Case I, we adopt the displacement and stress field of a semi- 
infinite elastic body B under a spherically-distributed pressure applied on a circular surface 
5c with radius a, i.e.. Love's solution based on the potential method [ ]. We then restrict 
attention to a spherical subbody >Bsph(i?) of radius R, slightly modified to accommodate the flat 
surface 5c (see Figure 20c). The surface points of interest to this analysis are give by (z, r) = 
(2i?sin^(6>/2), i?sin(6>)), with unit normal n = (- cos(6>), sin(6>)), for 9 G (arcsin(a/i?); tt]. 
Next, we replace the part of the semi-infinite body outside iBsph(i?) by the equilibrium surface 
traction q acting on 5/3sph(i?) ~ ^c, and we respect the displacement reference frame employed 
in Love's solution. 

Love's solution for a spherically-ditributed pressure applied on a circular surface with ra- 
dius a depends on the Newtonian potential V and the Boussinesq's logarithmic potential x 
[._]. Specifically, the displacements in cylindrical coordinates are 
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and the stress components of interest are 

(Jrr{z,r) 



1 

2^ 



2v—-(l-2v)^-z^ 
dz dr'^ dr'^ 






dV_ 9V 
dz dz'^ 
z d'^V 



27r drdz 
In these equations, the Newtonian potential and its derivatives are 
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and the derivatives of the Boussinesq's logarithmic potential are 
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with n given by the positive root of -^-^ + — = 1. 

Since the analytical evaluation of the surface traction q{&) and the normal component of 
the surface displacement Un{&) does not result in any further simplification of the equations 
presented above, we will only proceed to compare these results with those of the exact solution 
by means of numerical experimentation. We present this comparative study next in turn. 

Accuracy assessment 

The accuracy of the two approximate problems presented in this Appendix is assessed by com- 
paring their predictions with the exact solution recently presented by Zhupanska [ ]. Specifi- 
cally, an elastic incompressible sphere that accommodates a flat circular area of ajR = 0.233 
is employed. Since it is assumed that a <^ R, this example corresponds to a non-small contact 
area and therefore it allows for exploring the limits of validity of the presented formulation. 
Likewise, the selection of an incompressible material corresponds to the scenario with the 
largest displacements. 

The exact solution reported by Zhupanska [21] uses a general solution for the axisymmet- 
ric problem of an elastic sphere and a dual series equation approach to reduce the original 
boundary-value problem to a Fredholm integral equation of the second kind, which is then 
solved numerically. The analytical solution is achieved without replacing the contact area of 
the sphere by an elastic half-space and it allows for determining the contact pressure and the 
displacement of the traction-free surface. The results also show that the spherically-distributed 
contact pressure predicted by Hertz theory remains accurate for relatively large contact areas 
(i.e., the relative error in p)^ is smaller than 6% for a/R < 0.40, provided linear elasticity 
remains valid). 

Figure 21 shows the exact solution and the approximate solutions that correspond to the 
approximate problems referred to as Case I and Case II. The results are is very good agreement 
and also reveal Case I as being equally effective, and more efficient, than Case II. 
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Figure 21: Normal displacement of surface points for an elastic sphere under the action of 
contact stresses and supported at its center, for a/R = 0.233 and u = 0.50. The exact solution 
(symbols) is taken from Zhupanska [ ]. The approximate solutions correspond to approxi- 
mate problems referred to as Case I (black solid curve) and Case II (grey solid curve). The 
dashed curve corresponds to the assumption of traditional contact mechanics theories. 
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